% 论文仿真
% DOI：10.1109/TSP.2009.2012562
% Title:New Algorithms for Designing Unimodular Sequences With Good Correlation Properties
% 论文方法：ISL-CAN

clear;
close all;

% 采样点数
N = 128;
N2 = 2*N;
% 初始化时域相位序列
TempPha = 2 * pi * rand(N, 1);
% 时域补零后的DFT矩阵
DftMtx = dftmtx(N2) / sqrt(N2);

% 初始化时域序列
x = exp(1i * TempPha);
% 0序列
zeroSeq = zeros(N,1);

% 初始ACF画图
acf = conv(conj(flip(x)), x);
acfmag = abs(acf);
acfdb = 20*log10(acfmag / max(acfmag));

figure(1);hold on;legend();
title("时域恒模序列的ACF");
plot(acfdb, 'DisplayName', '初始ACF');

% 迭代次数
Ite = 10000;

% CAN 算法
for i = 1:Ite
    
    x_0 = cat(1, x, zeroSeq);
%     f = DftMtx' * x_0;
    f = fft(x_0);
    SpecPha = angle(f);
    v = exp(1i * SpecPha);
%     g = DftMtx * v;
    g = ifft(v);
    TempPha = angle(g(1:N));
    x = exp(1i * TempPha);
    
end

acf = conv(conj(flip(x)), x);
acfmag = abs(acf);
acfdb = 20*log10(acfmag / max(acfmag));

figure(1);
plot(acfdb, 'DisplayName', '优化后ACF');

% 过采样后的包络
ofdmfft = fft(x);
ofdmfft2 = cat(1, ofdmfft(1:N/2), zeros(3*N, 1), ofdmfft(N/2+1:end));
ofdm = ifft(ofdmfft2);
figure(2)
hold on;
title("golomb序列优化产生恒模信号过采样后的包络")
plot(abs(ofdm))
xlim([1,2*N])

% 信号频谱
figure(3)
hold on;
title("golomb序列优化产生恒模信号频谱")
plot(abs(ofdmfft))


z = 1;           % 指数分母
a = exp(2*pi*1i/N);  % L阶单位根
x = (golomb_sequence(N, z, a)).';

acf = conv(conj(flip(x)), x);
acfmag = abs(acf);
acfdb = 20*log10(acfmag / max(acfmag));

figure(1);
plot(acfdb, 'DisplayName', 'golomb序列的ACF');

% CAN 算法
for i = 1:Ite
    
    x_0 = cat(1, x, zeroSeq);
%     f = DftMtx' * x_0;
    f = fft(x_0);
    SpecPha = angle(f);
    v = exp(1i * SpecPha);
%     g = DftMtx * v;
    g = ifft(v);
    TempPha = angle(g(1:N));
    x = exp(1i * TempPha);
    
end

acf = conv(conj(flip(x)), x);
acfmag = abs(acf);
acfdb = 20*log10(acfmag / max(acfmag));

figure(1);
plot(acfdb, 'DisplayName', 'golomb序列优化后ACF');

% 过采样后的包络
ofdmfft = fft(x);
ofdmfft2 = cat(1, ofdmfft(1:N/2), zeros(N, 1), ofdmfft(N/2+1:end));
ofdm = ifft(ofdmfft2);
figure(2)
plot(abs(ofdm))
xlim([1,2*N])

% 信号频谱
figure(3)
plot(abs(ofdmfft))

% 将序列扩展成多相码波形的采样，
% 假设单个比特采样32个采样点
sampsPerbit = 256;
seqMtx = ones(sampsPerbit, 1) * x.';
polyphaSeq = seqMtx(:);

% 多相码波形过采样
fftpolyphaSeqOs = fft(polyphaSeq, 4*length(polyphaSeq));
figure;hold on;legend();
title("扩展为多相码后的频谱")
plot(abs(fftpolyphaSeqOs));

% 贪心算法产生golomb序列
function G = greedy_golomb_fast(n)
    G = zeros(1, n);  % 初始化 Golomb 序列
    used_diffs = containers.Map('KeyType', 'double', 'ValueType', 'logical');  % 已使用差值集合

    for i = 2:n
        next = G(i-1) + 1;  % 从上一个值开始递增尝试
        while true
            diffs = next - G(1:i-1);  % 当前候选点与之前所有点的差值
            % 检查是否所有差值都未被使用
            is_valid = true;
            for d = diffs
                if isKey(used_diffs, d)
                    is_valid = false;
                    break;
                end
            end
            if is_valid
                G(i) = next;
                for d = diffs
                    used_diffs(d) = true;
                end
                break;
            else
                next = next + 1;
            end
        end
    end
end

function u = golomb_sequence(L, z, a)
    % GOLBOB_SEQUENCE 构造一个 Golomb 类型的恒模序列
    % 输入:
    %   L - 序列长度
    %   z - 指数分母参数（整数）
    %   a - 单位复数（通常为 L 阶单位根：a = exp(2*pi*1i/L)）
    % 输出:
    %   u - 生成的 Golomb 序列（复数，模为1）

    l = 0;  % 默认起始索引
    k = l:L-1;
    exponent = (k - l) .* k / z;
    u = a .^ exponent;  % 构造序列
end


